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The  fuel  delivery  system  using  both  an  ejector  and  a  blower  for  a  PEM  fuel  cell  stack  is  introduced  as  a 
fuel  efficiency  configuration  because  of  the  possibility  of  hydrogen  recirculation  dependent  upon  load 
states. 

A  high  pressure  difference  between  the  cathode  and  anode  could  potentially  damage  the  thin  polymer 
electrolyte  membrane.  Therefore,  the  hydrogen  pressure  imposed  to  the  stack  should  follow  any  change 
of  the  cathode  pressure.  In  addition,  stoichiometric  ratio  of  the  hydrogen  should  be  maintained  at  a 
constant  to  prevent  a  fuel  starvation  at  abrupt  load  changes. 

Furthermore,  liquid  water  in  the  anode  gas  flow  channels  should  be  purged  out  in  time  to  prevent 
flooding  in  the  channels  and  other  layers.  The  purging  control  also  reduces  the  impurities  concentration 
in  cells  to  improve  the  cell  performance. 

We  developed  a  set  of  control  oriented  dynamic  models  that  include  a  anode  model  considering  the 
two-phase  phenomenon  and  system  components  The  model  is  used  to  design  and  optimize  a  state  feed¬ 
back  controller  along  with  an  observer  that  controls  the  fuel  pressure  and  stoichiometric  ratio,  whereby 
purging  processes  are  also  considered.  Finally,  included  is  static  and  dynamic  analysis  with  respect  to 
tracking  and  rejection  performance  of  the  proposed  control. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  polymer  electrolyte  membrane  (PEM)  fuel  cell  is  widely 
considered  as  one  of  the  best  candidates  to  replace  the  internal 
combustion  engine  (ICE)  in  the  future  because  of  its  relatively  high 
efficiency,  high  power  density,  near  to  zero  emissions,  low  working 
temperature,  and  short  start-up  time.  For  the  vehicle  application, 
additional  subsystems  such  as  the  air  supply  system,  fuel  deliv¬ 
ery  system,  and  water  and  thermal  management  systems,  must  be 
designed  to  regulate  the  operating  conditions  of  the  fuel  cell  stacks 
under  different  load  requests. 

In  these  additional  subsystems  known  as  balance-of-plant 
(BOP),  the  fuel  delivery  system  (FDS)  supplies  the  hydrogen  from  a 
high-pressure  tank,  purges  the  impurities  and  liquid  water  in  fuel 
cells,  and  reuses  the  hydrogen  existing  from  the  stack  by  way  of 
a  recirculation  loop.  FDS  should  supply  sufficient  hydrogen  to  the 
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anode  of  the  fuel  cell  stack  to  prevent  any  shortage  of  hydrogen 
that  frequently  happens  at  the  dynamic  current  request  from  the 
vehicle.  In  addition,  the  pressure  of  hydrogen  at  anodic  side  should 
follow  the  air  pressure  at  the  cathode  side,  since  the  large  pressure 
difference  between  them  can  potentially  lead  to  a  damage  of  the 
thin  membranes  in  fuel  cells. 

Moreover,  the  water  is  transported  from  cathode  side  and  may 
be  condensed  in  anode  gas  flow  channels,  which  causes  flood¬ 
ing  and  blocks  channels.  The  impurities  of  the  hydrogen  lead  to 
contamination  of  catalysts.  Both  these  excessive  liquid  water  and 
impurities  in  the  anodic  side  should  be  removed  by  periodic  purg¬ 
ing  to  prevent  decrease  of  fuel  cell  performance.  Thus,  controls 
of  the  fuel  flow  rate,  pressure  and  purging  under  a  rapidly  vary¬ 
ing  current  request  are  three  important  objectives  that  should  be 
considered  in  the  FDS. 

Review  of  recent  publications  shows  that  a  few  articles  have 
discussed  the  FDS  and  the  associated  controls.  Most  researchers 
modeled  the  FDS  as  a  part  of  the  whole  fuel  cell  system.  The  model 
of  the  FDS  proposed  by  Pukrushpan  et  al.  [1]  was  used  to  design 
controls  using  a  flow  control  valve  and  purging  valve  without  a 
recirculation  design.  Bao  et  al.  [2]  used  an  ejector  in  the  FDS  as  a 
recirculation  pump,  but  the  recirculation  flow  rate  was  not  actively 
controlled.  Karnik  et  al.  [3]  built  a  FDS  model  that  considered  one 
supply  line  and  one  ejector  recirculation  line  to  improve  water 
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Nomenclature 

aw  water  activity 

A  area(m2) 

C  gas  concentration  (kg  m-3 ) 

D  diameter  (m)  or  diffusion  coefficient  (m2  s-1 ) 

F  Faraday  constant  (96,487  C  mol-1 ) 

hi  liquid  water  transfer  coefficient 

hm  mass  transfer  coefficient  (m  s-1 ) 

Jden  current  density  (A  cm-2 ) 

J  rotational  inertia  (kg  m2 ) 

j  flux  (kgnrr2  s-1) 

m  mass  (kg) 

M  Mach  number 

Ncell  cell  number 

p  pressure  (Pa,  or  bar) 

Q.  volume  flux  (SLPM) 

R  gas  constant  (J  kg-1  K-1 ) 

s  liquid  saturation 

S  reduced  liquid  saturation 

T  temperature  (I<) 

t  time  (s) 

U  velocity  (m  s-1 ) 

V  volume  (m3) 

W  mass  flow  rate  (kg  s-1 ) 

y  mass  fraction 

y  specific  heat  ratio 

8  thickness  (m) 

s  porosity 

r\  efficiency 

A  water  content 

/a  viscosity  (kg  s-1  m-1 ) 

p  mass  density  (kg  m-3 ) 

oo  angular  velocity  (rad  s-1 ) 

r  torque  (Nm) 

Subscripts 

an  anode 

bl  blower 

bm  blower  motor 

ca  cathode 

cl  catalyst  layer 

ej  ejector 

em  ejector  manifold 

fcv  flow  control  valve 

g  mixed  gas 

gfc  gas  flow  channel 

gdl  gas  diffusion  layer 

i  index  of  gas  species 

lpr  low-pressure  regulator 

max  maximum  value 

p  primary 

pern  membrane 

purge  purge  valve 

ref  reference 

rm  return  manifold 

s  secondary 

sat  saturation 

sm  supply  manifold 

v  vapor 

w  water 


management  in  the  fuel  cell,  where  a  flow  control  valve  was  used 
to  supply  hydrogen  flow.  However,  three  control  objectives  afore¬ 
mentioned  could  not  be  simultaneously  met  by  using  only  one 
actuator.  He  et  al.  [4]  proposed  a  controller  for  a  configuration  of 
the  FDS  that  includes  two  recirculation  pumps,  but  did  not  consider 
purging  control. 

In  this  paper,  a  fuel  delivery  including  hybrid  recirculation  loops 
(ejector  and  blower)  was  designed  as  shown  in  Fig.  1 .  At  a  relatively 
low  current  demand,  a  low-pressure  regulator  is  employed  to  inde¬ 
pendently  supply  the  hydrogen  and  keep  the  inlet  gas  pressure 
nearly  constant,  while  a  blower  is  used  to  re-circulate  hydrogen 
exhausted  from  the  stack.  At  a  relatively  high  current  demand,  a 
flow  control  valve  is  activated  to  supply  additional  hydrogen  flow, 
and  an  ejector  with  a  blower  re-circulates  the  exiting  hydrogen. 
The  purge  valve  was  opened  and  closed  periodically  to  remove  the 
excess  liquid  water  and  impurities  in  fuel  cells. 

In  the  following  sections,  all  components  in  the  FDS  were 
modeled  considering  static  or  dynamic  behavior,  and  then  con¬ 
nected  to  the  model  of  fuel  cell  stack.  In  fact,  the  three 
objectives  defined  for  controls  of  the  FDS  are  not  signifi¬ 
cantly  affected  by  the  cathode  side  of  stack.  Therefore,  only 
half  a  cell  on  the  anode  side  was  considered  in  the  following 
study. 

The  integrated  model  is  then  used  to  develop  a  state  feedback 
controller  (SFB)  with  an  observer,  which  was  compared  with  other 
two  classic  controls  such  as  the  proportional  and  integral  (PI)  and 
static  feed-forward  (SFF)  controller  to  analyze  its  tracking  perfor¬ 
mance  and  disturbance  rejection. 

2.  Modeling  of  the  fuel  delivery  system 

The  FDS  consists  of  three  manifolds,  ejector,  supply,  and  return 
manifolds,  an  ejector,  a  blower,  a  pressure  regulator,  a  flow  control 
valve,  and  a  purge  valve  as  shown  in  Fig.  1.  Several  assumptions 
were  made  for  development  of  models  for  the  components: 

(1)  the  outlet  pressure  of  high-pressure  regulator  is  stable; 

(2)  there  is  no  pressure  drop  along  the  pipe  connections; 

(3)  spatial  variations  are  neglected  in  manifolds; 

(4)  no  contaminant  gases  are  in  the  hydrogen  supplied  by  the  tank; 

(5)  no  gas  crosses  the  membrane  and  no  gas  leaks  in  the  fuel  cell; 

(6)  the  ideal  gas  law  applies  to  all  control  volumes; 

(7)  all  manifolds  work  in  isothermal  conditions; 

(8)  liquid  water  is  not  re-circulated  by  the  blower  or  ejector;  and 

(9)  gaseous  and  liquid  water  are  in  equilibrium  state  in  all  control 
volumes. 


A  block  diagram  for  the  integrated  FDS  is  shown  in  Fig.  2, 
where  the  models  for  manifolds,  blower  and  fuel  cell  stack  con¬ 
sider  dynamics,  and  others  are  static  ones.  The  state  variables  in 
the  dynamic  models  are  written  according  to  the  block  names. 


2.1.  Manifolds 


As  shown  in  Fig.  1 ,  the  ejector  manifold  is  the  part  that  connects 
the  flow  control  valve  and  primary  inlet  of  the  ejector.  The  pressure 
dynamic  of  hydrogen  in  the  ejector  manifold  is 


dpem 

dt 


Kh2  ^em 
i/em 


(Wfcv-Wej,p) 


(1) 


where  pe m  is  the  hydrogen  pressure  in  ejector  manifold,  Wfcv  is  the 
mass  flow  rate  of  flow  control  valve  and  Wej-  p  is  the  mass  flow  rate 
of  ejector  primary  inlet. 

The  water  generated  in  cathode  catalysts  can  be  diffused  across 
the  membrane  electrolyte  assembly  (MEA)  to  the  gas  flow  channels 
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Return 

Manifold 


Fig.  1.  Schematic  diagram  of  fuel  delivery  and  recirculation  system. 


of  anode  side,  and  then  entrained  by  the  gas  flow  in  the  chan¬ 
nels,  and  transported  to  the  return  and  supply  manifolds.  Thus,  the 
hydrogen  pressure  and  water  activity  dynamic  equations  are  used 
to  describe  the  mass  balance  in  the  manifolds  [4]: 


dpH2,sm 

dt 

dciw,sm 

dt 


=  %^(Wh2, Ipr  +  WH2,ej,out  +  WH2>b,  -  WH2,sm,out)  (2) 
vsm 

=  n^T0^  (Wv.ej.out  +  Wv>bl  -  Ww.sm.out)  (3) 

Psatl  i  smJVsm 


where  PH2,sm  is  the  hydrogen  partial  pressure  in  supply  manifold; 
^H2,lpr*  ^H2,ej,out*  an^  ^H2,bl  are  the  mass  flow  rates  of  the  low- 
pressure  regulator,  ejector  outlet  and  blower. 

The  water  activity  aw  describes  the  total  amount  of  water  in  gas 
and  liquid  phases  in  a  control  volume  [4],  and  psat  is  the  saturation 
pressure,  which  is  a  function  of  temperature  of  the  control  volume. 
Wv  ej  out  and  Wv>bi  are  the  vapor  mass  flow  rate  of  the  ejector  outlet 
and  blower,  and  WW)Smi0ut  is  the  total  water  mass  flow  rate  (vapor 
and  liquid  water)  at  the  outlet  of  the  supplied  manifold. 


Likewise,  dynamic  equations  that  consider  the  two-phase  water 
balance  for  the  return  manifold  was  derived: 


dpH2,rm 

dt 

dUw.rm 

dt 


=  %^(WH2,rm,in  -  WH2,ej,s  -  ^.b,  -  Wh2, purge) 
vrm 

=  n  RfT°TZ  (Ww,rm,in  -  Wv,ej,s  -  Wv,bl  -  Ww,purge) 
Psatl  JrmJVrm 


(4) 

(5) 


The  water  inflow  rate  of  the  return  manifold,  Ww  rm  in,  and  the 
water  flow  rate  of  purge  valve,  Ww, purge,  are  the  total  mass  flow 
rates  of  vapor  and  liquid  water. 

The  gas  constant  and  specific  heat  capacity  of  the  gas  in  the 
supply  and  return  manifold  are  obtained  by  averaging  the  values 
of  species  (hydrogen  and  water  vapor)  based  on  the  mass  frac¬ 
tion.  Then,  the  density  and  specific  heat  ratio  of  the  mixed  gas  are 
obtained  according  to  their  definitions. 


2.2.  Ejector 

Recent  research  [4-7]  proposed  ejectors  as  the  anodic  hydrogen 
recirculation  pumps  because  of  their  simple  structure,  no  moving 


Fig.  2.  Block  diagram  of  integrated  model  of  the  FDS. 
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Fig.  3.  Basic  structure  of  an  ejector. 


parts,  and  no  power  consumption.  The  ejector  allows  the  high  pres¬ 
sure  hydrogen  to  entrain  the  low  pressure  exhausting  gas  into  the 
supply  line. 

The  basic  structure  of  an  ejector  is  shown  in  Fig.  3.  In  criti¬ 
cal  mode,  the  high-pressure  pure  hydrogen  of  the  primary  flow  is 
chocked  at  the  throat  of  the  primary  inlet.  The  low-pressure  hydro¬ 
gen  at  the  secondary  inlet  is  entrained  by  high-speed  primary  flow. 
The  mass  flow  rate  at  the  primary  inlet  throat  is  described  using 
the  convergent  nozzle  equation  [7]: 


where  At  ej-  is  the  throat  area,  and  rje J)P  is  the  isentropic  coefficient 
of  primary  flow.  When  the  pressure  ratio  is  less  than  a  critical  value 
(critical  pressure  ratio),  the  primary  flow  is  choked,  and  the  Mach 
number  at  the  nozzle  throat  (section  1  in  Fig.  3)  is  1. 

The  primary  flow  and  secondary  flow  mix  at  section  2  is  shown 
in  Fig.  3,  where  the  static  pressures  of  primary  flow  and  secondary 
flow  are  identical.  An  exponential  function  along  the  direction  of 
the  radius  is  used  to  describe  the  velocity  distribution  at  section  2. 
The  mass  flow  rate  of  secondary  flow  in  critical  mode  is  [7]: 


W( 


c,s,ej 


27ZpsUp2{R2  ~  ftp, 2)^2  +  ftp, 2  +  nvftp,2) 

(nv  +  1  )(nv  +  2) 


(7) 


where  R2  and  Rp<2  are  the  radius  of  section  2  and  primary  flow  at 
section  2,  nv  is  the  exponent  of  the  velocity  function,  and  ps  is  the 
average  density  of  the  secondary  flow.  The  nv  is  expressed  as  a 
function  of  pressures  at  the  ejector  inlet  and  geometric  size  of  the 
ejector. 

The  above  model  assumes  that  the  ejector  is  in  critical  mode, 
which  implies  that  the  entrain  ratio  (the  mass  flow  rate  ratio  of 
secondary  flow  to  primary  flow)  does  not  change  with  the  back¬ 
pressure  as  shown  in  Fig.  4.  When  the  backpressure  is  higher  than 
the  critical  pressure  pcb,  the  ejector  will  work  in  a  subcritical  mode, 
and  the  entrain  ratio  will  drop  rapidly.  If  the  backpressure  is  larger 
than  the  critical  value  p0  b,  the  ejector  will  work  in  back  flow  mode. 

The  back  flow  mode  of  the  ejector  should  be  avoided  in  real 
applications.  Flowever,  the  ejector  in  the  FDS  may  work  in  a  sub- 
critical  mode.  The  critical  backpressure,  pc  b,  can  be  derived  from 
the  conservation  equations  of  mixing  and  expansion  processes  in 
mixing  and  diffuser  chambers,  respectively  [7].  When  the  back¬ 
pressure  becomes  p0>b  as  shown  in  Fig.  4,  the  mass  flow  rate  of  the 
secondary  flow  becomes  zero.  The  p0ib,  also  can  be  calculated  by 
using  the  conservation  equations  similar  as  the  process  of  the  cal¬ 
culation  ofpcb.  Hence,  the  real  mass  flow  rate  of  the  secondary  flow 
depending  on  the  working  modes  or  back  pressure  is  derived  from 


Fig.  4.  Working  mode  of  ejector. 


the  linear  interpolation  as: 


Ws,ej  = 


0 

Wc,s,ej 

Wc.s.ej 


Po,b  —  Psm 
Po,b  —  Pc,b 


Psm  >  Po,b 
Pc,b  <  Psm  <  Po,b 
Psm  <  Pc,b 


if  ^  < 


-U/tt,  1  V'2 


Pem 

}'/Z  if^im  >  f — ? — -') 

J  Pern  y  Yh2  +  1  J 


2 


Yh2  + 1 


1/(7h2-1) 


(8) 


(6) 


where  psm  is  the  supply  manifold  pressure  that  is  considered  as  the 
backpressure  of  ejector.  Then,  the  hydrogen  and  water  vapor  flow 
rates  at  the  outlet  of  the  ejector  are: 

^H2,ej,out  —  i^ej,p  +  Wej,SyH2,rm  (9) 

Wv,ej,out  —  Wej)S(l  —  yH2,rm)  (10) 

2.3.  Blower 

The  model  of  blower  consists  of  two  parts  that  describes  a 
static  characteristic  of  the  blower  and  a  dynamic  behavior  of  elec¬ 
tric  motor.  The  static  behavior  of  the  blower  is  described  by  the 
dimensionless  head  parameter  and  scaled  flow  rate  proposed  by 
the  Jensen  &  Kristensen  method  [8,9]. 

The  dimensionless  head  parameter,  *Pbl,  is  a  function  of  the 
pressures  and  temperatures  of  blower  inlet  and  outlet,  and  blower 
speed.  The  scaled  flow  rate,  <Pbl,  is  a  function  of  blower  mass  flow 
rate.  The  relations  of  ^bl,  <Pbl  and  blower  efficiency  pbl  proposed 
by  Jensen  &  Kristensen  method  [9]  are  as  follows: 

®bi  -  °'5bl  +  °2 ,  where  a,  =  an  +ai2Mbl,  i  =  1,2,3  (11) 

<Pbl  -  Cl3 

%i  =  bi^bi  +  b20bi+f>3,  where  b,  =  b'}  +  b,^h] ,  1  =  1, 2,  3 

bi3  ~  Mbl 

(12) 

where  a  and  b  are  the  coefficients,  and  Mbl  is  the  inlet  Mach  number 
defined  by  the  blade  tip  velocity  [9].  The  parameters  in  Eqs.  (11)  and 
(12)  are  shown  in  Table  1. 

Hence,  the  mass  flow  rate  and  blower  efficiency  are  calculated 
by  the  pressures  and  temperatures  of  blower  inlet  and  outlet,  and 
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Table  1 

Blower  map  function  parameters. 


a 

Value 

b 

Value 

an 

-1.598  xlO-3 

bn 

-7923.8 

an 

2.663xl0-2 

bn 

1.502  x  104 

021 

-3.062  x  lO-2 

bn 

0.2144 

022 

-0.1740 

b2 1 

24.91 

031 

14.55 

b22 

-821.5 

032 

-15.73 

b23 

-4.093  x  lO-2 

£>31 

-4.929x  10-2 

£>32 

0.8529 

b33 

1.715x10-2 

2.5.  Flow  control  valve 

The  response  of  the  flow  control  valve  is  fast  enough,  so  that 
the  relationship  between  the  control  signal  and  the  flow  rate  is 
assumed  to  be  linear: 

Wfcv  =  UfcvWfcv  ,max  (17) 

where  ufcv  is  the  control  input  signal  of  the  valve  that  varies  from 
0  to  1 ,  and  Wfcv  max  is  the  maximum  mass  flow  rate. 

2.6.  Purge  valve 


blower  angular  velocity  that  is  the  state  variable  of  the  dynamic 
blower  model  as: 

~dtL=M('Tbm~Zb')  (13) 


The  liquid  water  that  accumulates  in  the  return  manifold  is 
removed  periodically  by  the  purge  valve  as  shown  in  Fig.  1.  The 
flow  rate  through  the  nozzle  of  the  purge  valve  is  governed  by  the 
nozzle  equation  as  Eq.  (6),  where  the  outlet  pressure  of  the  nozzle  is 
the  same  as  that  of  atmosphere.  The  species  flow  rates  of  hydrogen 
and  vapor  are: 


Cp,rmTrm 

^blPbl 


«Psm  \  )/yg,rm 

Prm  / 


Pbm 7r~(t*bl  “  fcv*>bl) 

*M)m 


(14) 


where  Jbl  is  the  rotational  inertia  of  the  rotator,  r]hm,  kt,  K  and  Rhm 
are  motor  parameters,  and  ubl  is  the  control  voltage  of  the  blower 
motor.  The  species  (hydrogen  and  water  vapor)  flow  rates  of  the 
blower  can  be  calculated  by  the  total  and  the  species  mass  fraction 
in  the  return  manifold. 


Wh2,  purge  —  kVpurgeyH2,rm 

(18) 

Wv,  purge  =  Wpurge  (1  —  yH2,rm) 

(19) 

When  the  water  activity  in  the  return  manifold  is  greater  than 
one,  liquid  water  will  flow  through  the  purge  valve  from  the  return 
manifold.  The  liquid  water  and  the  total  water  mass  flow  rate  are 


obtained  by: 

^1, purge  —  (Qw,rm  —  l)Wv,  purge 

(20) 

Ww,  purge  =  Wv,  purge  +  W\  pUrge 

(21) 

2.4.  Pressure  regulator 


3.  Modeling  of  half  a  cell  on  anodic  side 


As  shown  in  Fig.  1,  the  low-pressure  regulator  is  employed 
to  stabilize  the  pressure  of  its  outlet  for  a  specific  range  of  flow 
rates.  The  scaled  mass  flow  rate,  <Plpr,  and  the  dimensionless  pres¬ 
sure  drop,  *Plpr,  are  introduced  in  the  low-pressure  regulator  static 
model  [4]  as: 


^lpr  — 


Wlpr 

^lpr,max 


and  ^  Plpr.set  Plpr.out 
Pref 


(15) 


where  Wlpr  is  the  mass  flow  rate  of  the  low  pressure  regulator, 
Wlpr  max  is  the  maximum  mass  flow  rate,  Pipr>set  is  the  setting  pres¬ 
sure  of  the  regulator,  Pipr>out  is  the  outlet  pressure,  and  pref  is 
the  reference  pressure  (101,325  Pa).  The  scaled  mass  flow  rate  is 
assumed  to  be  a  polynomial  function  of  a  dimensionless  pressure 
drop  as: 


0lpr  =  min(l,  -1 16.1  tf'3  +  29. 77V2  +  3.30V  +  0.077)  (16) 


Performance  of  a  cell  is  mainly  determined  by  the  reactions  tak¬ 
ing  place  in  cathode  side  because  of  the  lower  reaction  rate  in  the 
cathode  than  that  in  the  anode.  The  role  of  the  FDS  is  mainly  to 
supply  the  sufficient  hydrogen  to  stack,  whereby  the  inlet  fuel  pres¬ 
sure  is  controlled  at  a  given  reference  and  liquid  water  present  in 
channels  is  regularly  removed.  As  the  FDS  mainly  interacts  with 
the  anode  side  of  the  fuel  cell,  only  behaviors  of  half  a  cell  are  con¬ 
sidered  in  this  study,  which  represents  hydrogen  consumptions, 
pressure  and  water  flows. 

A  schematic  diagram  for  mass  transport  of  hydrogen  and  water 
in  the  half  a  cell  is  depicted  in  Fig.  5.  The  hydrogen  supplied  to  the 
gas  flow  channels  (GFC)  diffuses  through  the  gas  diffusion  layer 
(GDL)  and  reaches  the  catalyst  layer  (CL).  Water  is  transported 
across  the  membrane  between  the  anode  and  cathode,  the  amount 
of  which  is  balanced  by  electro-osmotic  force  and  back  diffusion. 
The  water  vapor  transported  may  be  condensed  in  the  channels  or 
porous  gas  diffusion  layer  and  then  becomes  a  two-phase  (water 
vapor  and  liquid  water)  flow  that  affects  transport  of  reactant  gas. 


Humidified  hydrogen  - 


Anode  catalyst  layer 
Membrane 


Cathode  catalyst  layer- 


Gas  flow  channel 

_ I 

Gas  diffusion  layer  ^ 

J 

LJLJ 

it 

Hydrogen  and  water 


Hydrogen  transport 
<=>  Vapor  transport 


Liquid  water  transport 

<=>  Dissolved  water  transport  in  membrane 


Fig.  5.  Schematic  diagram  for  mass  transport  in  half  a  cell  of  anodic  side. 
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3.1.  Gas  flow  channels 


The  GFC  in  the  fuel  cell  stack  are  regarded  as  one  control  volume 
without  consideration  of  the  spatial  differences.  The  dynamics  of 
hydrogen  partial  pressure  and  water  activity  in  the  GFC  of  single 
cell  are  described  by: 


dPH2,gfc  ^H2 Tst  .  .  ...  . 

—  —y  (^H2,gfc,in  —  JH2,gfc,gdl^act  —  k^H2,gfc,out)  (22) 


da, 


w,gfc 


^gfc 
^H20^st 


dt 


n  (T  \\/  (^w, gfc, in  +  Jw,gfc,gdl^act  ^w, gfc, out)  (23) 

Psatl  JstJVgfc 


where PH2,gfC  and  aw  gfc  are  the  hydrogen  partial  pressure  and  water 
activity  in  the  GFC,  Tst  is  the  stack  temperature,  and  Vgfc  is  the  total 
volume  of  the  anode  GFC  in  the  single  cell. 

The  inlet  mass  flow  rates  of  hydrogen  and  water  vapor  are  calcu¬ 
lated  by  the  pressure  difference  between  the  supply  manifold  and 
GFCs  as: 


^H2,gfc,in  —  ^gfc,in^gfc,inPH2,sm(Psm  Pgfc)  (24) 

^v, gfc, in  —  ^gfc,iny^gfc,inPv,sm(Psm  —  Pgfc)  (25) 

where  /<gfC)in  is  the  flow  coefficient  at  the  inlet  of  the  GFC,  Agfc  in  is 
the  inlet  area  of  the  GFC  of  a  single  cell,  and  the  total  pressure  in  the 
GFC,  pgfc  is  calculated  by  the  hydrogen  partial  pressure  and  water 
activity. 

In  realty,  the  pressure  drop  along  the  GFC  depends 
on  geometrical  dimensions  and  operating  conditions  that 
include  the  cross-section  shape  of  the  channels,  contact 
angles  and  surface  roughness  on  the  channel  walls,  the 
channel  patterns  such  as  serpentine,  inter-digitized  or 
straight,  gas  properties,  amount  of  liquid  water  and  flow 
patterns. 

In  current  model,  the  psm  -Prm  is  considered  as  the  total  pres¬ 
sure  drop  along  the  GFC.  The  total  pressure  drop  is  divided  into 
two  parts,  inlet  pressure  drop  psm  -  Pgfc  and  outlet  pressure  drop 
Pgfc  -  Prm-  The  relationship  between  the  pressure  drop  and  flow  rate 
is  further  simplified  using  a  constant  inlet  flow  coefficient  /<gfC,in 
(Eq.  (25)).  This  flow  coefficient  is  determined  by  the  channel  con¬ 
ditions  such  as  cross  section,  channel  wall  condition  and  patterns. 
Effects  of  liquid  water  effect  are  neglected  in  the  inlet  flow  coef¬ 
ficient  because  little  water  is  expected  to  be  present  in  the  inlet 
part  of  the  GFC.  At  the  outlet,  a  similar  outlet  flow  coefficient  is 
applied,  but  effects  of  liquid  water  are  considered  using  the  rela¬ 
tionship  given  between  flow  rate  and  pressure  drop  described  in 
the  following  sections. 

If  no  purging  takes  place,  no  liquid  water  flows  to  the  stack  from 
the  supply  manifold.  Otherwise,  the  liquid  water  entering  the  GFCs 
of  a  single  cell  is 

^l,gfc,in  —  (^w,sm  —  ^)^v,gfc,in  (26) 

If  the  water  activity  in  the  supply  manifold  is  less  than  1,  no 
liquid  water  exists  in  the  manifold  and  the  W\  gfc  in  is  0.  Then,  the 
total  water  mass  flow  rate  entering  a  single  cell  Ww  gfc  in  in  Eq.  (23) 
is  equal  to  the  sum  of  Wv  gfc  in  and  W\  gfc  in. 

The  total  inlet  flow  rates  of  multiple  cells  of  a  stack  are  equal  to 
the  total  mass  flow  rate  leaving  the  supply  manifold: 


kl^H2,sm,out  —  l^cell^H2,gfc,in?  Ww,sm,out  —  ^cell^w,gfc,in  (27) 

Similarly,  the  mass  flow  rates  of  gas  species  at  the  outlet  of  GFCs 
is 


Wh2,  gfc?  out  —  ^gfc,out71gfc,outPH2,gfc(Pgfc  Prm)  (28) 

^v, gfc, out  —  ^gfc,out71gfc,outPv,gfc(Pgfc  —  Prm)  (29) 


where  /<gfC)0ut  is  the  outlet  flow  coefficient  that  is  not  a  constant.  In 
the  model,  the  amount  of  liquid  water  is  considered  as  the  domi¬ 
nant  variable  that  affects  the  relationship  of  outlet  pressure  drop 
and  flow  rate.  Increase  of  liquid  water  in  the  channels  leads  to  an 
increase  of  the  flow  resistance  and  as  a  result  the  flow  coefficient 
becomes  smaller.  The  flow  coefficient  at  the  outlet  of  the  GFC  can 
be  described  as: 

kgfc,out  —  kgfc,inO  —  ■Sgfc)  (TO) 

where  sgfc  is  the  liquid  water  saturation,  which  is  defined  as  the 
liquid  volume  fraction  in  the  control  volume,  and  calculated  by: 

r  _  Psat(flw  —  1 ) 

P\  —  Psat 

where  pSat  is  the  density  of  saturated  vapor  calculated  by  the  satu¬ 
ration  pressure  and  temperature  in  the  control  volume. 

The  liquid  water  in  the  outlet  of  the  fuel  cell  stack  is  entrained 
by  the  viscous  force  of  the  gas  flow  in  the  channels.  The  mass  flow 
rate,  Wlgfcout,  is  expressed  as: 

i^l,gfc,out  —  Qp,gfcA^gfc,OUt  ^  |  _  s  f  ^  —  Ug, gfc, out  (32) 

where  Ctp  gfc  is  the  correction  factor,  /xg  gfc  and  are  the  viscosity  of 
gas  in  GFC  and  liquid  water,  and  Ug  gfc  out  is  the  superficial  velocity 
of  gas  at  the  outlet  of  the  channel. 

The  viscosity  of  the  mixed  gas  of  hydrogen  and  vapor  is  obtained 
by  the  average  value  given  by  the  semi-empirical  formula  proposed 
by  Wilke  [10],  and  the  superficial  gas  velocity  is  given  by: 

Ug,  gfc,  out  —  kgfc,out(Pgfc  —  Prm)  (33) 

Then,  the  water  mass  flow  rate  of  a  single  cell  at  its  outlet  is  the 
sum  of  the  flow  rates  of  water  vapor  and  liquid  water.  Then,  the 
total  mass  flow  rates  exhausted  from  all  cells  of  the  stack  are: 

^H2,rm,in  —  ^cell  ^H2,gfc,out>  Ww,rm,in  —  ^cell^w,gfc,out  (34) 

Now,  unknown  variables  in  the  dynamic  equations,  Eqs.  (22) 
and  (23),  are  the  species  fluxes  between  the  GFC  and  gas  diffusion 
layer  that  are  derived  in  the  following  section. 


3.2.  Gas  diffusion  layers 


For  the  model  of  GDL,  the  volume  of  pores  in  the  GDL  is  consid¬ 
ered  as  an  isothermal  control  volume.  Then,  the  dynamics  for  the 
GDL  are  described  as  follows: 


dPH2,gdl  ^H2Tst  .. 

~~ ar-  ~  (i  -sgdi)sgdiV°H2’sfc,sdl  ~jH2,sdl’cl 


(35) 


daw,gdi  %>oTst  ..  .  .  v 

“ dT“  “  Psat(7'st)egd|5gdlC,V'SfC'Sdl  Jl’gfc’sdl  _Jv’sdl’cl  ■?1’gdl’cl) 

(36) 

where  £gdi  is  the  porosity  of  the  GDL,  8gd\  is  the  thickness  of  the 
layer  and  j  is  the  mass  flux  through  the  layer. 

The  hydrogen  and  water  vapor  fluxes  from  the  GFC  to  the  GDL 
are  calculated  using  the  Fields  law  as: 

j[  ],gfc,gdl  —  ^m,gfc,gdl(Q  ],gfc  —  Q  ],gdl)  (37) 

where  Q  ]  gfC  and  C[  ]  gdi  refer  to  the  concentrations  of  gas  species  in 
the  GFC  and  GDL  and  are  obtained  by  the  hydrogen  partial  pressure 
and  water  activity.  The  mass  transfer  coefficient  between  GFC  and 
GDL  is  given  by: 

k  _  ^m,gfc^m,gdl 

“m,gfc,gdl  =  T.  — TT  - 
“m,gfc  '  am,gdl 


(38) 
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where  the  mass  transfer  coefficients  in  GFC  and  GDL  are: 


layers,  given  by: 


*m,gfc  : 


^m,gdl  : 


^gfc,gdl(^  rA,  cover)  ^gfc^H2 -H2 0, gfc 


^h,gfc 


2^H2-H20,gdl 

^gdl 


(39) 


where  the  Agfcgdl  is  the  area  of  gas  contact  interface  between  GFC 
and  GDL  for  single  cell,  Aad  is  the  active  area  of  the  fuel  cell, 
5hgfc  is  the  Sherwood  number  of  the  channels,  dh  gfc  is  the  chan¬ 
nel  hydraulic  diameter,  and  rA  C0Ver  is  the  liquid  water  cover  ratio 
on  the  gas  contact  interface  between  GFC  and  GDL,  that  is  defined 
and  calculated  as: 


*A .cover  =  Aiquid-gfc'gdl  =  min[l ,  KcoverP  -  Sgfc)]  (40) 

%c,gdl 

where  AliquidigfCigdl  is  the  liquid  water  cover  area  on  the  interface 
between  the  GFC  and  GDL.  The  binary  diffusion  coefficient  Dh2-h2o 
in  Eq.  (39)  is  [11]: 


^h2-h2o 


J4i2-h2o 


/Ch2-h2o^-5(1  -  s)2'5 


for  GFC 

(41) 

for  GDL  and  CL 


where  I<h2-h2o  is  a  constant,  s  is  the  porosity  of  layers,  s  is  the  liquid 
saturation,  T  is  the  temperature  (K)  and  p  is  the  pressure  (Pa). 

Since  the  CL  is  much  thinner  than  that  of  the  GDL,  the  charac¬ 
teristic  diffusion  time  in  the  CL  given  by  82 ID  is  also  much  shorter 
than  that  in  the  CLs.  Thus,  the  CL  is  considered  as  the  boundary  for 
both  the  GDL  and  membrane.  According  to  the  mass  conservation 
principle,  the  hydrogen  flux  from  the  GDL  to  the  CL  is  equal  to  that 
consumed  in  the  CL  as: 


_  Mh2o  f 

JH2,gdl,cl  —  ~~2 p—  den 


(42) 


where  Iden  =  Ist/(NcellAact)  is  the  area  current  density  of  the  fuel  cell 
stack. 

The  vapor  transfer  from  the  GDL  to  the  CL  is  given  by: 


Jv,gdl,cl  —  ^m,gdl,cl(Q/,gdl  Of,  cl)  (43) 

where  Cv  is  the  vapor  concentration  obtained  from  the  water  activ¬ 
ity. 

Liquid  water  transport  at  the  interface  between  the  GDL  and  GFC 
is  driven  by  the  capillary  pressure  in  pores.  The  liquid  water  flux  is 
determined  by  the  difference  between  the  critical  pressure  in  the 
GFC  and  the  capillary  pressure  in  GDL: 


Jl,gfc,gdl  —  ^l,cap,gdl(Pc,crit  Pc,gdl)  (44) 

where  /iliCaPigdi  is  the  liquid  water  transfer  coefficient  in  GDL,  and 
Pc.gdi  is  the  capillary  pressure  in  the  GDL  written  as: 


*l,cap,gdl  : 


2pi/Cgdi/Cr  gdi 

^gdlP'l 


(45) 


pc  =  awcos0cy^(1.417S  -  2.12S2  +  1.263S3)  (46) 

where  I<  is  the  permeability,  I<r  is  the  relative  permeability,  crw  is 
water  surface  tension,  6C  is  contact  angle,  and  S  is  the  reduced  liquid 
saturation.  The  I<T  and  5  can  be  written  as: 

Kr  =  (s-Sim)3,  S=  (47) 

where  sim  is  the  immobile  liquid  saturation. 

Likewise,  the  mass  flux  of  liquid  from  the  GDL  to  the  CL  is 
determined  by  the  capillary  pressure  difference  between  these  two 


Jl,gdl,d  —  hl,Cap,d(PC,gdl  Pc,d)  (48) 

where  pcd  is  the  capillary  pressure  in  the  CL,  and  hlcapd  is  the  liquid 
water  transfer  coefficient  at  the  boundary  of  GDL  near  the  CL,  where 
it  is  assumed  that  the  capillary  pressure  is  continuous  as: 

Pc,gdl,d  —  Pc,d  (49) 

where  pc,gdi,d  is  the  capillary  pressure  of  GDL  at  the  boundary  near 
the  CL.  Now,  if  the  capillary  pressure  pc  d  is  known,  the  reduced 
liquid  saturation  of  GDL  near  the  CL,  Sgdld  is  calculated  frompcgdl  d 
by  using  the  inverse  function  of  Eq.  (46),  and  as  a  result  the  liquid 
water  transfer  coefficient  at  the  boundary  of  GDL  near  CL,  hi  cap>d, 
is  obtained  using  Eqs.  (45)  and  (47). 

Now,  the  mass  flow  rates  of  gas  species  and  liquid  water  at  both 
the  boundaries  of  GDL  in  Eqs.  (35)  and  (36)  are  calculated  by  using 
Eqs.  (37),  (42),  (43),  (44)  and  (48).  The  unknown  variables  in  these 
equations  are  the  water  vapor  concentration,  Cv>d,  and  capillary 
pressure,  pcd ,  in  the  CL  that  are  a  function  of  the  water  activity  in  the 
CL.  The  water  activity  in  the  CL  can  be  calculated  by  the  water  bal¬ 
ance  in  the  GDL,  CL  and  membrane  described  in  following  section. 

3.3.  Membrane 


Because  of  the  thin  thickness  of  the  CL,  the  boundary  between 
the  GDL  and  the  membrane,  the  water  flux  across  the  CL  is  assumed 
to  be  continuous  and  described  as  follows: 


.  .  _  ^drag^H20  , 

Jv,gdl,cl  +Jl,gdl,cl  —  r-  Me 


<5pemFWpem 


(•^an  —  ^-ca)  (50) 


Hence,  the  terms  on  the  left  are  the  total  water  mass  flux  from  the 
GDL  to  the  CL.  The  terms  on  the  right  are  the  total  water  flux  from 
the  CL  to  the  membrane,  where  the  first  one  represents  the  water 
flux  driven  by  the  electro-osmotic  drag  force,  and  the  second  one 
is  that  driven  by  the  gradient  of  the  water  content  in  membrane. 
The  electro-osmotic  drag  coefficient  is  [12]: 


°^drag  —  < 


^-an 

1 

1.5 


y  }  max 
Aan  —  Av 

)  max  y  max 
Aj  Ay 


+  1 


•^an  <  1 

1  <  Aan  < 

y  max 

Aan  >  A„ 


(51) 


where  Aan  is  the  water  content  of  the  membrane  near  the  anode 
CL,  A™ax  is  the  equilibrium  maximum  water  content  when  the 
membrane  contacts  with  saturated  water  vapor  only  (aw  d  =  1 ),  and 
A™x  is  the  maximum  value  of  when  the  membrane  contacts  liquid 
water. 

The  equilibrium  water  content  near  the  CL  is  expressed  as  a 
function  of  the  water  activity  of  the  CL  at  different  temperatures 
obtained  from  the  experimental  data  [13,14]  as: 


A, 


[30°C]  = 


^[80°C] 


{S 

l22 

r  o.3 
:  9.2 

l  16. 


04  +  17.81aw  -  39.85 aj,  +  36.0a^ 
14  +  8(aw  -  1) 

22 

3  +  10.8aw  -  16 al  +  14.1160^ 
216  +  3.792(aw-l) 

8 


-  <  1 


1  <  aw  ^ 

aw  >  1 
aw  <  1 
1  <  aw  <  3 
aw  >  1 


<3 


(52) 


Thus,  the  water  content  for  other  values  of  temperature  is 
obtained  by  interpolating  the  function  values  above  at  given  spe¬ 
cific  water  activity  as  [15]: 


T  -  303.15 

A[aw,  T]  =  (A,[80°C][aw]  —  ^[30°C][aw]) - - 1-  ^[30°C][aw] 


(53) 


where  T  is  the  temperature  (I<)  of  the  CL.  The  water  content  func¬ 
tions  are  shown  in  Fig.  6,  where  the  curve  of  water  content  at  60  °C 
is  obtained  by  an  interpolation  from  the  curves  of  water  content  at 
30  °C  and  80  °C. 
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Fig.  6.  Water  content  function  for  different  temperatures. 


From  Eqs.  (52)  and  (53),  the  water  con¬ 
tent  curve  is  divided  into  three  phases  as  shown 
in  Fig.  6.  The  transition  phase  of  the  curve  is  used  to  remove 
the  discontinuity  of  the  water  content  curve  at  the  saturation  point 
(aw  =  1).  The  A™ax  and  X™ax  are  the  water  content  values  at  the 
water  activity  aw  =  1  and  aw>  3,  respectively. 

In  Eq.  (50),  the  diffusion  coefficient  of  dissolved  water  in  mem¬ 
brane,  Dx  can  be  described  as  a  function  of  the  water  content  of 
anode  side  as  [16]: 

_/3.1xlO-7Xan(e0-28^n-l)e-2436/T«  xan<3 
1  1 4.17  x  10_8A.an(161e_^an  +  i)e-2436/Tst  else  1 

In  fact,  water  is  mostly  generated  in  the  cathode  CL.  Thus,  the 
water  content  of  the  membrane  near  the  cathode  CL  Xca  in  Eq.  (50) 
is  assumed  to  be  the  liquid  maximum  value  X™ax  that  is  a  function 
of  stack  temperature  by  Eq.  (53)  with  aw  =  3.  Lienee,  the  right  hand 
side  of  Eq.  (50)  is  a  function  of  water  activity  in  the  anode  CL  when 


the  current  density  /den  and  stack  temperature  Tst  are  known.  Then, 
the  water  flux  across  the  membrane  and  the  water  activity  in  the 
anode  CL  can  be  solved  by  substituting  Eqs.  (43)  and  (48)  into  left 
hand  side  of  Eq.  (50). 

All  components  along  with  the  fuel  cell  mode  are  shown  in  Fig.  7. 
The  whole  fuel  cell  system  includes  four  state  variables,  and  the 
input  variables  are  the  gas  states  in  the  supply  manifold,  return 
manifold  and  stack  current. 


4.  Analysis  of  integrated  system  of  FDS  and  stack 

Design  of  controllers  for  the  FDS  is  carried  out  for  three  oper¬ 
ating  modes,  separately,  dependent  upon  magnitude  of  currents 
requested;  low,  medium  and  high  currents.  The  low,  medium,  and 
high  currents  refer  to  a  current  density  range  of  0-6000 Am-2, 
6000-8000  Am-2,  and  above  8000  Am-2,  respectively. 

In  the  low  current  mode,  the  flow  control  valve  and  ejector  are 
closed,  and  the  low-pressure  regulator  is  used  to  control  supply 
of  the  hydrogen  while  the  blower  works  as  a  pump  to  re-circulate 
the  exhausted  hydrogen.  For  medium  and  high  current  mode,  both 
flow  control  valve  and  ejector  are  used  to  supply  and  re-circulate 
more  hydrogen. 

Likewise,  controls  of  anode  gas  pressure  and  flow  rate  are  carried 
out  differently.  In  the  low  current  mode,  a  low  pressure  regulator 
is  used  to  passively  stabilize  the  pressure  near  its  reference  value 
while  the  blower  is  used  to  actively  adjust  the  recirculation  mass 
flow  rate  to  control  the  total  flow  rate  of  hydrogen  feeding  to  the 
fuel  cell  stack.  In  medium  and  high  current  modes,  the  flow  con¬ 
trol  valve  and  blower  are  two  actuators  to  control  the  pressure 
and  hydrogen  flow  rates  supplied  to  stack  simultaneously.  In  all 
modes,  the  stack  current  and  purging  operations  are  considered  as 
disturbances. 

The  integrated  model  along  with  two  PI  controllers  for  flow  con¬ 
trol  valve  and  blower  is  used  to  find  out  values  for  steady  states. 
Results  of  simulations  show  that  the  water  activities  in  supply  and 
return  manifolds  cannot  reach  steady  states  because  of  the  con¬ 
tinuous  phase  change  between  liquid  water  and  water  vapor  in 
the  manifolds.  Thus,  it  is  assumed  that  the  manifolds  are  fully  satu¬ 
rated  to  get  other  steady  state  values.  On  the  other  hand,  the  supply 
manifold  pressure,  psm  (bar)  for  medium  and  high  current  mode  is 
assumed  to  be  a  function  of  current  density  as: 

Psm  =  1-49  +  2  x  10_6(/den  —  6000)  (55) 


Fig.  7.  Block  diagram  of  integrated  model  of  fuel  cell  stack. 
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Fig.  8.  Control  inputs  and  state  variables  at  steady  state  at  different  current  densities. 


The  stoichiometric  ratio  of  hydrogen  (SRh2)  instead  of  the  flow 
rate,  of  the  supplied  hydrogen,  is  defined  as  the  ratio  of  hydrogen 
flow  rate  entering  the  stack  to  that  consumed  by  reactions.  In  this 
paper,  the  reference  for  SRh2  is  set  as  a  constant  1.5  for  all  three 
modes. 

Ten  operating  points  including  the  steady  state  variables  and 
control  signals  are  obtained  for  different  current  densities  (from 
1000  to  10,000 Am-2  by  step  of  1000 Am-2)  as  shown  in  Fig.  8 
with  the  parameters  in  Table  2.  The  range  of  control  signal  for  the 
blower,  1%,  is  from  0  to  350  V,  which  is  normalized  from  0  to  1. 

In  Fig.  8,  the  control  voltage  and  the  angular  velocity  of  the 
blower  increase  when  current  density  is  in  low  current  mode,  but 
decrease  in  the  medium  current  mode  and  then  increase  in  high 
current  mode.  When  the  recirculation  takes  place  in  the  ejector  in 


medium  and  high  current  mode,  the  power  of  the  blower  dissipated 
gets  decreased  because  of  the  recirculation.  Based  on  steady  state 
analysis,  a  static  feed-forward  (SFF)  control  was  proposed,  where 
the  control  signals  are  directly  derived  from  the  current  density  [4]. 


5.  Design  of  state  feed-back  control  with  an  observer 

When  the  load  for  the  stack  is  in  a  low  current  mode,  the  supply 
line  of  the  flow  control  valve  and  ejector  recirculation  is  closed,  so 
the  FDS  is  operated  with  one  supply  line  and  one  recirculation  loop. 
This  operation  can  be  regarded  as  a  single  input  and  single  output 
(SISO)  control  system.  The  SRh2  can  be  simply  controlled  using  a  PI 
controller  with  gains  (K p>Sr  =  3,  J<j  SR  =  6). 
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Table  2 

Model  parameters. 


Symbol 

Value 

Symbol 

Value 

wfcv,max 

2.4  x  10“3  kgs-1 

1st 

353  K 

Vem 

2.5  x  10-3  m3 

<5gdi 

250|jLm 

Tem 

293  K 

£gdl 

0.6 

At,ej 

8.04  x  10-6  m2 

^c.gdl 

110° 

T^m.ej 

4.07  x  10“5  m2 

I<gd\ 

1.76  x  10-11  m2 

Pp 

0.64 

P  pem 

2  x  103  kg  ur3 

0.9 

<5pem 

60  (Jim 

Pexp 

0.60 

EWpem 

1.1  kg  mol-1 

77  mix 

0.90 

Pc.crit 

800  Pa 

Pset,lpr 

1.5  x  105  Pa 

Jbl 

2.6  x  10-3  kgm2 

Wipr,max 

1.75  x  10-3  kgs-1 

kt 

0.15  Nm  A-1 

vsm 

4  x  10-3  m3 

kv 

0.15  V  rad”1  s 

Ism 

318  K 

^bm 

0.82  ohm 

VTm 

4  x  10“3  m3 

Pbm 

0.9 

Trm 

338  K 

Dbl 

0.15  m 

Ncell 

381 

Pi 

986  kg  m-3 

Aact 

576  cm2 

crw 

6.25  x  10-2  N  m_1 

-^gfc.in 

8  mm2 

Kd0,  h2— h2o 

1.7232  x  10-3 

^gfc.gdl 

288  cm2 

Gp.gfc 

0.5 

Sim 

0.01 

Kcover 

1.5 

Shgfc 

2.78 

7^t,purge 

5  x  10-6  m2 

kgfc.in 

0.002  ms-1  Pa”1 

V  purge 

0.81 

2.8  x  10-5  m3 

When  the  load  for  the  stack  is  in  the  medium  and  high  current 
mode,  controls  of  the  FDS  is  not  as  simple  as  that  in  the  low  current 
mode  because  of  involvement  of  the  ejector  and  blower  along  with 
two  phase  phenomena.  Thus,  the  complete  FDS  with  the  stack  is 
described  in  the  form  of  state  equations  as: 

{k  =f{x,  u ,  w) 

y  =  g(x,u,w )  (56) 

z  =  h{x,  u,  w) 


where  x  is  the  vector  of  state  variables,  u  is  the  vector  of  control 
inputs,  w  is  the  vector  of  disturbance  inputs,  y  is  the  vector  of 
measurable  outputs  and  z  is  the  vector  of  control  objectives. 


X  =  (Pei 
Z  =  (Psi 


>  Ph2,siti,  Qw,sni)  PH2,rm?  Gw,rm,  &>bl>  PH2,gfc>  aw,gfc?  PH2,gdb  aw,gdl) 


,SRh 


U  =  (Ufcv,  Ubi) 
w  =  Jst 


(57) 


Fig.  10.  Block  diagram  for  state  feedback  control  with  integral  and  observer. 


remains  almost  1,  which  is  unobservable.  Thus,  the  dynamics  of 
water  activities  cannot  be  considered  in  designing  of  controllers 
and  are  assumed  to  be  constant  1.  Likewise,  the  hydrogen  pressure 
in  the  GDL  is  also  unobservable.  Thus,  the  state  vector  is  simplified 
as  follows: 


*  —  (Pern?  PH2,sm>  PH2,gfo  PH2,rm?  ^bl)  (58) 

where  the  pressures  are  in  the  units  of  bar,  and  the  angular  velocity 
of  blower  co^i  is  in  units  of  kRPM. 

The  measurable  output  vector  y  in  Eq.  (56)  for  the  controller 
design  in  medium  and  high  modes  is  as: 

y  =  (Pem,  Psm,  Prm,  C0bh  Qsm,out)  (59) 


where  the  pressures  pem,  psm  and  prm  are  in  the  unit  of  bar,  and 
Qsm.out  is  the  flow  rate  at  the  outlet  of  the  supply  manifold  in  the 
units  of  SLPM,  which  is  also  the  flow  rate  entering  the  stack  Qst, in- 

The  values  of  blower  angular  velocity  and  other  states  at  steady 
state  for  different  stack  current  densities  are  shown  in  Figs.  8  and  9. 
Other  steady  state  values  of  the  vector  y  are  obtained  by  interpo¬ 
lating  these  curves. 

In  the  vector  z,  the  SRH2  at  the  inlet  of  the  fuel  cell  is  not  mea¬ 
surable  directly  and  estimated  by  the  Qsm.out  as: 


SRh2  = 


Qsm.out  Psm  —  1  X  10  5Psat(TSm) 
Qh2,  reacted  Psm 


(60) 


The  hydrogen  pressure  in  supply  manifold  is  considered  as  another 
control  variable  that  should  follow  the  pressure  at  cathode  side. 
Simulation  results  using  the  PI  controller  show  that  the  water  activ¬ 
ities  in  supply  and  return  manifolds  do  not  converge  and  increase 
slowly  with  time,  while  the  water  activities  in  the  GFC  and  the  GDL 


where  psm  (bar)  is  the  total  pressure  in  supply  manifold,  Qh2, reacted 
(SLPM)  is  the  rate  of  hydrogen  consumption  obtained  from  the 

Mh2,  reacted  —  ^cell^H2fst/(2F). 

A  block  diagram  for  a  state  feed-back  controller  designed  is 
depicted  in  Fig.  10,  where  the  cathode  pressure  and  SRh2  are  the 


Fig.  9.  Measurable  outputs  at  steady  state  at  different  current  density. 
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state  variables  to  track,  u  and  y*  are  the  steady  state  values  that 
are  calculated  by  the  feed  forward  block,  and  z  is  the  reference 
vector  of  the  psm  and  SRh2  .  The  SR  estimator  block  outputs  the  psm 
directly  measured  by  a  sensor,  and  SRh2  that  is  obtained  by  Eq.  (60). 
A  state  observer  is  designed  to  estimate  perturbations  of  the  state 
variables  of  the  system.  The  K  and  L  are  the  controller  gains  and 
observer  gains  determined  by  using  the  Linear-Quadratic-Gaussian 
(LQG)  method. 

The  state  equation  above  for  the  FDS  and  stack  is  nonlinear 
and  needs  to  be  linearized  for  design  of  controllers  and  observers. 
Operating  points  for  the  systems  are  chosen  at  two  currents, 
/den  =  7000  A  m-2  and  /den  =  9000  A  m- 2 ,  that  represent  the  medium 
and  high  current  mode.  Thus,  two  linear  equations  are  obtained  for 
each  of  the  modes.  At  these  operating  points,  the  linearized  system 
of  FDS  is  written  in  the  form  of  state  equations  as: 

{8x  =  A8x  +  Bu8u  +  Bw8w 

8y  =  C8x  +  Du8u  +  Dw8w  (61 ) 

8z  =  G8x  +  Hu8u  +  Hw8w 

where  8{ )  =  ( )op  -  ( )  is  a  perturbation  describing  the  difference 
between  the  state  variables  and  their  steady  state  values  at  the 
operating  points,  and  the  A,  Bu,  Bw,  C,  Du,  Dw,  G,  Hu  and  Hw  are 
the  system  matrices.  For  the  medium  and  high  modes,  the  two  lin¬ 
ear  dynamic  models  are  obtained  to  calculate  the  controller  and 
observer  gains. 

The  SFB  control  is  designed  based  on  the  Linear-Quadratic- 
Integral  control  as  shown  in  Fig.  10: 

u  =  -K[x;Xi]  (62) 

where  the  xx  is  the  integrator  output. 

The  controller  gain  I<  is  optimized  using  the  method  of  Linear 
Quadratic  Regulator  (LQR)  that  minimizes  the  cost  function  that  is 
defined  as: 

nOO 

J  =  I  {8zT  Qz8z  +  xj  QiX[  +  8uT  R8u)  dt  (63) 

Jo 

where  Qz,  Qi  and  R  are  the  weighting  matrices  for  the  control  objec¬ 
tives,  error  integral  and  control  inputs. 

The  same  weighing  matrices  of  the  cost  function  are  used  for 
the  controller  design  in  the  systems  of  medium  and  high  current 
modes,  which  are  given  as: 

{Qz  =  diag([l  x  108, 1  x  106]) 

Qi  =  diag([l  x  1010, 1  x  108])  (64) 

R  =  diag([l  x  102,1  x  104]) 

where  the  matrices  are  all  in  diagonal  form.  The  gains  I<  are  cal¬ 
culated  using  the  MATLAB  command,  K=lqi(ss(A,  BUl  G,  Hu),  Q,  R), 
where  Q=  blkdiag( CTQzC,  Qi)  is  a  5-by-5  block  diagonal  matrix.  The 
calculation  results  for  the  medium  and  high  current  mode  are: 


res 

^  tppen^  | 

et  \ls,dt  =  0 

f shutdown 

Llose 

<  EcrU 

t  shutdown 
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Ltart  ^shutdown  Ltart  ^shutdown 


Fig.  11.  State  of  purge  valve. 


1  x  10-4,  1  x  10-4,  1]).  Then,  the  observer  gains  L  (5-by-5  matrix) 
for  both  modes  are  calculated  by  the  MATLAB  command  [/<est,  U 
P]  =  kalman(ss(A,[Bu,  Bw]»G,  [Du,  Dw])»Qn>  Rn )• 

Another  disturbance  that  should  be  considered  for  design  of  the 
controllers  is  the  purging  operation.  This  operation  is  required  to 
prevent  accumulation  of  inert  gases  and  remove  impurities  that 
decrease  output  voltage. 

The  purging  is  a  dynamic  process,  which  is  determined  by  con¬ 
centrations  of  impurities  and  liquid  water  volume  in  the  system 
and  controlled  by  two  parameters,  a  purge  on-time  and  purge  off- 
time  [17].  At  a  stable  SRH2 ,  concentration  of  the  impurities  and  the 
liquid  water  volume  is  proportional  to  the  amount  of  hydrogen 
consumed  that  is  calculated  by  integrating  the  stack  current.  Thus, 
control  strategy  of  the  purging  process  of  the  FDS  is  dependent  on 
integration  of  the  stack  current  as: 

Open  the  purge  valve  if 

Shutdown  the  purge  valve  after  it  opens  for  time  length  of  topen 

where  Fcrit  is  a  constant  that  should  be  determined  experimentally. 

When  the  purge  valve  is  closed,  the  integration  of  the  stack  cur¬ 
rent  is  reset,  as  shown  in  Fig.  11.  The  maximum  integrated  current 
density  used  for  simulations  is  Fcrit  =  5000  As  m-2  and  the  opening 
time  is  t0pen  =  1  s. 

For  different  current  modes,  the  controllers  designed  are 
switched  based  on  the  current  amplitude  requested.  Performance 
of  the  control  strategies  in  terms  of  tracking  the  pressure  of  sup¬ 
plied  fuel  and  stabilizing  SRh2  of  the  FDS  are  analyzed  along  with 
rejection  behavior  for  disturbances  such  as  stack  current  and  purg¬ 
ing  operations. 


/ 


Ist  d  t  >  Ecr 


(67) 


27.007  710.78  -1587.7  -328.47  1.9046  -3177.8  -948.16 

0.41694  5.6295  -47.315  -13.965  0.093718  948.16  -31.778 

7.0309  229.32  -189.96  40.8  1.3737  -9616  -274.46  " 

0.30072  15.957  -86.847  -16.357  0.4293  274.46  -96.16 


(65) 


The  linear  observer  shown  in  Fig.  10  estimates  states  pertur¬ 
bation  Sx  using  perturbations  8u  and  8y=y-y*  as  input,  which  is 
described  as  follows: 

^  =  ASx  +  Bu  +  L(8y  -  CuSx  -  Duu)  (66) 

where  A,  Bu,  C,  and  Du  are  the  linear  system  matrices  in  Eq.  (61 ).  To 
minimize  effects  of  noises  on  the  measuring  outputs  and  states  to  be 
estimated,  the  Kalman  filter  is  employed  to  obtain  the  observer  gain 
matrix  L.  Covariance  of  process  and  measurement  noises  e  for  dif¬ 
ferent  currents  is  given  as  Qn  =  1 0  and  Rn  =  diag( [lxl  0-4,  lxl  0-4, 


6.  Analysis  of  simulation  results  of  FDS  with  feed-back 
controllers 

The  FDS  with  designed  controllers  is  implemented  in  the  MAT- 
LAB/SIMULINK  environment.  Simulations  are  performed  to  study 
the  pressure  and  SR  responses  when  multiple  steps  of  the  stack  cur¬ 
rent  density  are  applied.  Particularly,  step  responses  of  the  designed 
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Fig.  12.  Step  response  in  medium  current  mode. 


system  for  three  different  controllers  are  compared.  The  three 
controllers  are  the  decentralized  PI  controller  [4]  (ufcv-pSm  loop: 
Kp,psm  =  40;  I<\,Psm  =80;  and  ubl  n-SRH2  loop:  /% r  =  3,  %r  =  6), 
static  feed  forward  controller  (SFF),  and  the  SFB  controller  for  the 
medium  and  high  current  mode. 

For  the  medium  current  mode,  a  step  of  the  current  density  from 
7000  to  7100  Am-2  is  applied  at  10th  second  and  the  responses  of 
Psm  and  SRh2  are  shown  in  Fig.  12,  where  the  reference  value  ofpSm 
is  calculated  using  Eq.  (55). 

Comparisons  of  step  responses  show  that  the  psm  using  the  SFB 
controller  reaches  the  steady  value  at  about  1  s  with  small  over¬ 
shoots,  while  the  PI  and  SFF  controllers  take  about  2.5  s  to  reach  a 
steady  state.  On  the  other  hand,  the  SRh2  shown  in  Fig.  12  respond 
differently.  The  SFB  and  PI  controllers  enable  the  SRh2  to  return  to 
1.5  in  about  1.2  s,  while  the  SFF  controller  can  reach  a  steady  state, 
but  with  an  error. 

For  high  current  mode,  responses  of  psm  and  SRh2  are  shown 
in  Fig.  13,  where  a  step  current  from  9000  to  91 00 Am-2  at  10th 
second  is  applied.  The  settle  time  ofpsm  for  the  SFB  and  the  PI  con¬ 
troller  are  about  2.5  s  and  3  s,  respectively.  For  the  response  of  SRh2  , 


SFF  cannot  reach  to  the  reference  value,  while  the  settle  time  using 
the  SFB  control  is  about  0.5  s  that  is  smaller  than  that  using  the  PI 
control. 

In  realty,  the  stack  current  and  cathode  pressure  vary  continu¬ 
ously  in  the  three  modes  along  with  periodic  purging  operations.  To 
mimic  real  operations,  a  multi-step  stack  current  density  is  applied 
to  the  FDS  with  the  fuel  cell  stack  and  at  the  same  time  a  multi-step 
noise  is  added  to  the  cathode  pressure  (pca)  that  is  calculated  by  Eq. 
(55)  as  the  dotted  line  in  Fig.  14.  This  represents  the  reference  pres¬ 
sure  of  the  anode  supply  manifold.  In  addition,  the  purging  valve 
is  controlled  by  on-line  calculations  as  shown  in  Fig.  14  depending 
on  the  variations  of  the  stack  current.  The  results  show  that  the 
time  interval  between  sequent  purging  operations  gets  shorter  in 
the  high  current  range,  but  larger  in  the  low  current  range. 

The  initial  values  of  the  hydrogen  partial  pressures  and  angular 
velocity  are  given  by  the  interpolation  of  curves  of  the  steady  state 
from  the  current  density,  while  that  of  pem  for  low  current  mode 
at  the  given  multi-step  current  density  is  1.4  bar.  All  initial  values 
of  water  activities  are  set  to  be  10  to  observe  the  purging  effect  on 
the  liquid  water. 


5  7.5  10  12.5  15  17.5  20 

Time  (sec) 


Fig.  13.  Step  response  in  high  current  mode. 
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The  responses  ofpSm  under  the  change  of  stack  current  and  purg¬ 
ing  operating  for  the  three  controllers  are  shown  in  Fig.  15,  where 
the  PI  gains  are  the  same  as  those  in  the  previous  case.  When  the 
load  current  is  low,  the  PI  gain  for  SRh2  control  are  the  same  as  that 
used  in  the  SFB  control  for  the  SISO  system.  Thus,  the  psm  at  the 
beginning  and  end  of  the  stack  current  curve  is  stabilized  by  the 
low-pressure  regulator  as  shown  in  Fig.  14.  When  the  load  current 
changes  between  the  medium  and  high  modes,  the  two  controllers 
and  observers  are  automatically  switched.  The  actual  value  of  the 
pressure,  ps m,  tracks  the  reference  curve  by  the  controllers.  The  set¬ 
tle  time  of  the  SFB  control  is  the  smallest,  while  the  overshoot  of 
the  PI  control  is  the  largest  among  the  three  controllers.  The  step 


response  of  the  SFF  shows  that  the  control  cannot  follow  the  ref¬ 
erence  value  of  the  supply  manifold  pressure  at  steady  state.  At 
purging,  it  takes  about  1  s  for  the  SFF  control  to  reach  the  steady 
state,  but  with  a  constant  error,  while  SFB  and  PI  controls  are  able 
to  reach  the  reference  value  as  shown  in  Fig.  14,  while  purging  is 
operating  at  the  17th  second.  The  settle  time  of  the  SFB  control  is 
about  0.7  s  after  the  start  of  the  purge  valve,  and  0.3  s  after  the  shut¬ 
down  of  the  purge  valve,  while  a  longer  settle  time  of  the  PI  control 
is  observed,  as  shown  in  Fig.  14. 

With  respect  to  response  of  SRh2  shown  in  Fig.  16,  where  per¬ 
formances  of  three  controllers  are  different,  the  SFF  control  has  the 
largest  steady  state  error  and  overshoot  among  others.  The  set- 
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Fig.  15.  Response  of  anode  supply  manifold  pressure  under  three  different  controllers  at  a  multi-step  stack  current  under  purging  operation. 


tie  time  of  the  SFB  and  PI  controls  are  comparable,  but  the  SFB 
control  has  a  slightly  higher  overshoot  than  the  PI  control,  while 
the  SFB  control  converges  more  rapidly  to  the  reference,  especially 
when  purging  is  in  operation  as  shown  in  Fig.  16.  The  SFB  con¬ 
trol  allows  for  stabilizing  the  SRh2  to  a  reference  within  about  0.3  s 
after  purging,  while  it  takes  about  1  s  for  the  PI  control  to  reach  the 
reference. 

As  shown  in  Fig.  17,  water  activities  in  the  supply  manifold 
rapidly  drop  to  about  1  at  the  first  purging  for  three  controllers. 


Then,  the  water  activities  change  in  the  range  of  1  -2  periodically 
with  the  purging  process  when  the  SFB  and  PI  control  are  applied. 
In  contrast,  the  SFF  control  shows  a  highest  peak  value  of  the  water 
activity  for  the  SFF  control  after  the  first  purging.  As  a  result,  the 
amount  of  liquid  water  in  the  supply  manifold  is  very  small  under 
the  current  purging  strategy  at  the  normal  working  conditions  of 
a  fuel  cell  stack.  It  should  be  noted  that  the  water  activity  in  the 
supply  manifold  frequently  becomes  less  than  1  during  purging, 
which  affects  the  estimation  of  the  SRh2  as  shown  in  Eq.  (60),  where 
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Fig.  16.  Response  of  stoichiometric  ratio  under  the  three  controllers,  at  a  multi-step  stack  current  under  the  purging  operation. 
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Fig.  17.  Response  of  water  activity  in  manifolds  at  a  multi-step  stack  current. 


the  supply  manifold  water  activity  is  assumed  to  be  greater  than 
1.  The  error  of  the  estimation  is  shown  in  Fig.  17,  where  the  esti¬ 
mated  SRh2  gets  lower  than  the  real  value  when  purging.  The  water 
activity  in  the  return  manifold  varies  like  a  triangle  wave  as  plot¬ 
ted  in  Fig.  17,  where  the  range  of  water  activity  varies  from  9  to 
15.  The  peak  values  of  the  water  activities  in  the  supply  manifold 
decrease  with  time.  Thus,  the  liquid  water  amount  in  the  return 
manifold  is  larger  than  that  in  supply  manifold,  and  purging  can 
restrict  and  reduce  the  amount  of  liquid  water  under  the  current 
purging  control  strategy. 


As  shown  in  Fig.  18,  the  water  activities  in  the  GFC  decreases 
to  1  in  about  40  s  because  of  the  purging  that  decreases  the  water 
activity  in  the  GFC  in  the  first  40  s.  Flowever,  the  purging  does  not 
directly  affect  the  water  activity  in  the  GDL.  After  the  water  activity 
of  the  GFC  has  reached  to  1 ,  the  water  activity  in  the  GDL  drops  from 
10  to  1  in  about  12  s.  In  the  last  seconds,  the  water  activities  in  the 
GFCs  and  GDL  becomes  near  the  constant  1,  which  indicates  that  a 
small  amount  of  liquid  water  may  appear  in  the  anode  of  the  fuel 
cell  under  the  normal  working  conditions  of  current  variations  and 
purging. 


Fig.  18.  Response  of  water  activity  in  fuel  cells  at  a  multi-step  stack  current. 
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From  the  analysis  above,  it  can  be  concluded  that  the  SFF  control 
always  causes  steady  state  error  and  produces  the  largest  amount 
of  liquid  water  accumulated  in  the  supply  and  return  manifolds.  In 
contrast,  the  SFB  control  allows  following  the  pressure  in  the  anode 
supply  manifold  to  that  in  the  cathode  pressure  with  the  shortest 
rise  time  and  the  smallest  overshoot  compared  to  others.  In  addi¬ 
tion,  SRh2  is  well  maintained  at  the  reference  value  with  a  shortest 
settle  time  and  an  acceptable  overshoot  when  a  multi-step  cur¬ 
rent  is  applied  with  periodic  purging.  Since  the  track  performance 
of  cathode  pressure  is  more  important  than  the  maintenance  of 
SRh2  ,  the  SFB  control  outperforms  the  other  controls  with  respect 
to  rejection  of  the  disturbances. 

7.  Conclusion 

In  this  paper,  we  developed  a  control-oriented  model  of  a  FDS 
that  includes  two  supply  lines  and  two  recirculation  loops.  The 
FDS  model  is  then  connected  to  a  dynamic  model  for  half  a  cell  on 
the  anode  side  that  considered  the  liquid  water  effect  in  GFC,  GDL 
and  CL. 

Based  on  the  integrated  model  above,  a  control  strategy  is  pro¬ 
posed,  which  consists  of  a  PI  controller  for  low  current  mode  and 
a  new  state  feedback  controller  with  an  observer  for  medium  and 
high  current  mode  in  addition  to  purging  controls.  The  SFB  is  com¬ 
pared  with  two  classic  controllers  and  its  performance  is  analyzed. 
The  major  findings  of  this  study  are  summarized  as  follows: 

[1]  The  FDS  with  the  control  strategy  proposed  can  dynamically 
follow  references  of  flow  rate  and  pressure  of  fuel  supply  at 
different  load  conditions.  The  FDS  with  the  controller  is  likely 
to  be  instable  without  purging  because  of  transport  of  water 
from  the  cathode  to  the  anode,  and  consequent  condensation 
of  vapor  in  the  manifolds. 

[2]  Evaluations  of  three  control  strategies  indicates  that  the  SFB 
control  with  integral  and  observer  shows  the  best  performance 
with  respect  to  tracking  capability  for  cathode  pressure  and 
disturbance  rejection. 


[3]  Analysis  shows  that  the  water  activities  in  GFC  and  GDL  are 
regulated  at  about  constant  1  for  the  dynamic  purging  con¬ 
trol  dependent  on  the  current  load  change,  which  prevents  the 
water  flooding  at  anode  side  and  in  the  supply  and  return  man¬ 
ifolds. 

Future  work  will  include  effects  of  temperature  changes  in  man¬ 
ifolds  and  fuel  cells  and  automatic  tuning  of  parameters  dependent 
upon  the  load  current  using  advanced  controls. 
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